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Abstract. An off-lattice, continuous space Kinetic Monte Carlo (KMC) algo- 
rithm is discussed and applied in the investigation of strained heteroepitaxial 
crystal growth. As a starting point, we study a simplifying (l + l)-dimensional 
situation with inter-atomic interactions given by simple pair-potentials. The 
model exhibits the appearance of strain-induced misfit dislocations at a char- 
acteristic film thickness. In our KMC simulations we observe a power law de- 
pendence of this critical thickness on the lattice misfit, which is in agreement 
with experimental results for semiconductor compounds. We furthermore in- 
vestigate the emergence of strain induced multilayer islands or Dots upon 
an adsorbate wetting layer in the so-called Stranski-Krastanow (SK) growth 
mode. At a characteristic kinetic film thickness, a transition from monolayer 
to multilayer island growth occurs. We discuss the microscopic causes of the 
SK-transition and its dependence on the model parameters, i.e. lattice misfit, 
growth rate, and substrate temperature. 



1. Introduction 

Atomistic models of crystal growth in Molecular Beam Epitaxy and similar tech- 
niques continue to attract considerable interest, see for example |21 El 0] for 
overviews, a very brief introduction is given in 0, this volume. Many interesting 
aspects of epitaxial growth can be discussed for the case of homoepitaxy where 
only one material is present. This includes the deposition of, for instance, a metal 
or a semiconductor compound on a matching substrate crystal. 

The term heteroepitaxial crystal growth refers to situations, where the de- 
posited adsorbate material differs from the substrate. The mismatch can be quite 
fundamental, as for instance in the growth of organic films upon metal substrates. 
Perhaps the most frequent and, certainly, the conceptually simplest situation is 
that where adsorbate and substrate would crystallize in the same type of lattice, 
but with slightly different lattice spacings. 
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The latter case is of particular relevance in the fabrication of semiconduc- 
tor heterostructures for potential technical applications, such as modern opto- 
electronic or storage devices. Prominent examples are the deposition of Ge on Si, 
InAs on GaAs, or CdTe on ZnSe. Frequently, the aim is to produce a smooth ad- 
sorbate film of a well defined thickness on a given substrate. However, the misfit 
induced strain can perturb the lattice structure and result in, e.g., dislocations 
which reduce the quality of the film drastically. 

On the other hand, mechanisms of strain relaxation can also be exploited 
in epitaxial growth. The most prominent example is perhaps the self-organized 
formation of multilayered islands in the so-called Stranski-Krastanov growth mode 
1 . The ideal process would yield dislocation-free three-dimensional objects of well 
defined size and shape which furthermore display spatial ordering. These so-called 
Quantum Dots can capture single or few electrons and may play the role of artificial 
atoms in a novel type of solid-state-lasers, for example. For a recent overview of the 
many interesting experimental and theoretical aspects of Quantum Dot physics, 
see [§]■ 

Besides the obvious technological relevance of heteroepitaxial growth, it is 
also highly interesting from a theoretical point of view. The self-assembly of Dots 
might be considered a prototype example of self-organized ordering processes. 
Heteroepitaxy remains furthermore a challenge in the design of growth models 
and the development of novel simulation techniques. 

The natural tool for mismatched heteroepitaxy appears to be Molecular Dy- 
namics (MD) jj] or quantum mechanical variations thereof, e.g. (HJ- This very 
attractive but computationally demanding approach is discussed in the general 
context of epitaxial growth in [5] (K. Albe, this volume), see ^U] f° r an example 
application in the context of dislocation formation. Whereas modern computers 
allow for the simulation of very large systems, the most serious limitation of the 
MD technique remains: The physical time intervals that can be targeted are gen- 
erally quite small, e.g. on the order of picoseconds. MBE relevant time scales of 
seconds or even minutes do not seem feasible currently, even when applying highly 
sophisticated acceleration techniques as in 

In simplifying lattice gas models, the particles can only be placed exactly at 
pre-defined sites. Nevertheless, it is possible to incorporate certain misfit effects 
into such models. In fact, some of the earliest KMC simulations were performed 
in the context of strained heteroepitaxial crystal growth El- So-called ball and 
spring models consider additional harmonic interactions, i.e. springs, between ad- 
sorbate or substrate atoms at neighboring lattice sites, see El and EI f° r 
more recent examples. In other approaches, continuum elasticity theory is applied 
to evaluate the contribution and influence of strain for a given configuration of the 
lattice El (and references therein). 

However, traditional lattice based models cannot cope with the emergence of 
defects or dislocations in a straightforward fashion. In the following we will discuss 
an off-lattice model and the corresponding Kinetic Monte Carlo (KMC) technique, 
which allows for continuous particle positions and in which strain results directly 
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from the atomic interactions. Schindler and Wolf have essentially formulated the 
model and simulation method in |171 118j , and some of the results presented here 
have been published previously in |191 121)1 12T] . For a detailed description of the 
concept and its applications see also |22j . 

The method incorporates the essential features of hetero-epitaxy in a natural 
way, and, at the same time, enables us to perform simulations over reasonable 
physical time scales. In contrast to earlier realizations of similar concepts, see e.g. 
|171 1181 1231 , we are also able to simulate the growth of rather thick films over 
a wide range of misfits. 

In section |5| we will define the model and describe the simulation technique. 
We briefly discuss the role of misfit dislocations as a mechanism for strain relax- 
ation and present our simulation results in section [3] Section 0] summarizes our 
investigations of the Stranski-Krastanov growth mode and in section[S]we conclude 
with a brief outlook on perspective investigations. 

2. Model and method 

In the following, we address quite fundamental and general aspects of heteroepi- 
taxy, rather than the properties of specific material systems. Hence, the primary 
goal is to develop a fairly simple model which still captures the essential features 
of a mismatched growth situation. 

Following several previous studies ^1 EI HJO I2H > we choose a pair potential 
ansatz to model interactions between the particles. To begin with, we consider 
a Lennard- Jones (LJ) system [2E]- A strong repulsive term at small distances 
competes with an attractive contribution in the potential 



The relative distance of particles i and j varies continuously with their position 
in space. By choice of the parameters U and a we can specify the different material 
properties in our model: interactions between two substrate or adsorbate atoms are 
given by the sets {U s ,a s } and {U a ,a a }, respectively. In principle, an independent 
set of parameters would define the interaction of substrate with adsorbate atoms. 
In order to keep the number of parameters small, we follow a standard approach 
and set U as = ^/U s U a , a as = {cs + &a) /2 for the inter-species potential. 

The potential energy of two isolated atoms that interact via (|2.1I) is minimal 
for rij = 2 1//6 cr. Accordingly, the lattice spacing in a Lennard- Jones crystal is 
proportional (and close) to a as well. The relative lattice misfit in our model is 
therefore directly controlled by choice of a s and a a : 



Essential features and several principled questions can already be studied in the 
simplifying framework of (l + l)-dimensional growth, where particles are deposited 
on a one-dimensional surface. For the simulations presented here we have prepared, 




(2.1) 



e = (a a - cr s ) /a s . 



(2.2) 
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Figure 1. Illustration of the Molecular Statics calculations in 
(1 + 1) dimensions. As the (black) adatom is virtually moved to 
the right, along the x-direction, the total energy of the system is 
locally minimized with respect to all other particle coordinates, 
including z for the moving adatom. The resulting potential energy 
surface displays an essentially periodic shape along the terraces 
and a pronounced Schwoebel barrier for moves across the edge. 

at least, 6 layers of substrate particles, each containing L particles. The system 
sizes vary between L = 200 and L = 800 in the following, depending on the 
computational costs of the considered problem. We assume periodic boundary 
conditions in the lateral direction and fix the particle positions in the bottom 
layer in order to stabilize the crystal. 

In our simulations of the growth kinetics, adsorbate particles are deposited 
at a constant rate Rd as measured in monolayers per second {ML/ s). Desorption 
events will be completely suppressed here as they would be very unlikely anyway. 
Only adsorbate particles at the surface are considered mobile in the sense that they 
perform diffusion hops and the like. However, bulk crystal adsorbate and substrate 
atoms will be due to relaxation processes which modify their spatial positions and 
relative distances without changing the topology. 

The most distinctive feature of the method is the treatment of thermally 
activated processes, for which we assume Arrhenius like rates, see ^ El E]- 
In lattice gas models, particles are moved from one site to another with a rate 
which usually depends on a small neighborhood only. In the off-lattice approach 
any event concerns the entire configuration and its rate depends on all particle 
coordinates, in principle. At a given time, the system resides in a local minimum 
of the total potential energy of the n particles: 

n n 

E *ot = E E u *i- ( 2 - 3 ) 

i=l .7=1+1 

A possible event k takes the system to one of the neighboring minima in config- 
uration space and for any such change the corresponding energy barrier and 
rate i?& has to be evaluated by means of a so-called Molecular Statics calculation, 
see [161 126] for examples and further references. 
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In a most general setting the identification of potential events, i.e. the relevant 
energy minima, can be difficult |26l 1271 125]. The (1 + l)-dimensional situation 
simplifies these calculations to a large extent. 

Figure n illustrates the procedure: The single particle on top of the terraced 
surface is virtually moved to the right, along the ^-coordinate in the illustration. 
For every fixed value of x, the total energy 12.. 'ill can be (locally) minimized with 
respect to the moving atom's z-coordinate and the positions of all other atoms 
including the substrate. In doing so, one can evaluate a potential energy surface 
(PES) as displayed in Fig. ^ Of course it is not necessary to scan the PES con- 
tinuously for the Arrhenius rates, as only the energies in the initial state and the 
transition state are relevant. The latter corresponds to a maximum in the figure, 
i.e. a saddle point in the energy as a function of all coordinates. In more complex 
situations like (2 + l)-dim. growth it is also possible to determine saddle points 
in a PES by iterative gradient based methods |22l 1261 1771 I28| . However, the iden- 
tification of the physical path from one minimum to the next can be much more 
involved in general, e.g. for concerted moves in (2 + 1) dimensions. 

Note that the (1 + l)-dim. system displays a strong Ehrlich-Schwoebel effect 
0E]: hops from a terrace are hindered by a barrier which relates to the very weakly 
bound transition state right at the edge. The effect is also present, but much weaker 
in a (2 + l)-dim. fcc-system, for instance. There, a descending particle can follow 
a more favorable path between neighboring particles at the edge. 

For large misfits, exchange processes at terrace or island edges may become 
more frequent than downward hopping diffusion. We have checked that for the 
values of e considered in the following, exchange diffusion can be neglected, in 
general. The precise barriers, however, depend in a subtle way on the misfit, the 
island size, and the properties of the actual interaction potential. For details we 
refer to [22] and forthcoming publications. 

Further simplifications may be used here. First, because the LJ-potential (|2.1(l 
is very weak for distances > 3<r, one can effectively cut-off the interactions and 
consider only particles in a corresponding neighborhood region when evaluating 
the PES of a particular event. The so-called frozen crystal approximation should 
be employed with care: in this scheme all coordinates but those of the moving 
atom are considered fixed, which significantly simplifies the saddle point search. 
Although the scheme appears to impose a drastic restriction, it has been reported 
by various authors that its main effect in the LJ-system is a uniform shift of all 
barriers by roughly 10% [IHl 113 CHI H2 • 

In our simulations, exchange diffusion and multiple jumps have been ne- 
glected. We consider only hops to the left or right neighboring minimum in the 
PES. Given a configuration of the system, we set up a catalogue of all these diffu- 
sion events and evaluate their rates 



R k = v Q exp 



E}; 



with 



v = 10 12 /s 



(2.4) 



k B T 



with Ek obtained as described above. 
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Figure 2 . Dislocations: a film grown in the model without down- 
hill tunneling JHj for a large misfit of e = +10%, section of a 
larger system with six layers of substrate. The lighter a particle is 
displayed, the smaller is its average distance to the nearest neigh- 
bors. Dislocations form right at the interface of adsorbate and 
substrate, as marked by the additional filled circles. Note that 
the elastic interaction with the film also affects the substrate. 

For simplicity, we choose a constant pre-factor or attempt frequency v Q for 
all diffusion processes, see for a brief discussion and references. Next, one of 
the diffusion or deposition processes is chosen with a probability proportional to 
its rate. In order to facilitate an efficient book-keeping and fast search, events are 
stored in a binary search tree, see [2| and references therein. After performing 
the particular process, the catalogue of rates has to be updated and re-evaluated. 
Obviously this task consumes a large portion of the computing power. Nevertheless, 
the rejection- free method is advantageous over a rejection-based simulation 0. 

Approximations and numerical inaccuracies, and even more so the use of the 
frozen crystal picture, result in deformations of the crystal and the accumulation 
of artificial strain energy. In order to avoid this effect, a relaxation of the system 
should be performed after each event, in the sense that the configuration is taken 
to the nearest local minimum of E to t- In practice we perform after each diffusion 
or deposition a local relaxation within a radius 3ct about the location of the event. 
The full global procedure is applied only after a certain number of events, which 
is decreased whenever the relaxation changes the barriers by more than, say, a few 
percent. 

In the following sections we discuss applications in the context of two different 
mechanisms of strain relaxation. Their peculiarities require certain extensions and 
modifications. However, the basic concepts and essential ingredients have been 
outlined above. 



3. Formation of dislocations 

If the misfit is relatively small in heteroepitaxial growth, one frequently observes 
the initial growth of a pseudomorphic film. Here the substrate determines the 



Heteroepitaxial growth 



7 



lateral lattice constant also in the adsorbate and strain energy accumulates in the 
growing film. On the other hand, for very thick adsorbate films far away from 
the substrate, one clearly expects to find a lattice structure and spacing as in 
the undisturbed bulk material. The ultimate relaxation of strain is through misfit 
dislocations, lattice defects which allow the adsorbate to assume its natural lattice 
structure, eventually. Here we address the question of how the transition from 
strained pseudomorphic to relaxed adsorbate growth occurs and how the lattice 
spacing evolves with the film thickness. 

In the frame of our model we choose the same potential depth for all types 
of interactions, which defines the energy scale. By identifying U a = U s = U as with 
the value U = 1.3l25eV, for instance, we obtain a barrier of 0.90eV for surface 
diffusion in the case of zero misfit. The substrate temperature was set to T = 
0.03U o /kB ~ 460if which is well below the melting temperature of a monoatomic 
Lennard- Jones crystal. All results presented in this section were obtained with the 
deposition rate = lML/s. 

Dislocations can be identified by constructing Voronoy polyhedra and search- 
ing for particles with a coordination number different from six in the bulk. Further 
information is obtained from a Burgers construction, see |191 \22\ for details. 

In our earlier investigations presented in |19| we considered positive and neg- 
ative values of the misfit parameter in the range —15% < e < +11%. For large 
absolute values of e the system displays the formation of misfit dislocations right at 
or very close to the substrate/adsorbate interface, cf. Fig. [21 The slightly oversim- 
plifying picture is that, initially, separated islands or mounds grow on the substrate 
and dislocations emerge where they meet. These are then overgrown by the de- 
posited material. Diffusion of dislocations by concerted moves of the surrounding 
particles are highly improbable in the LJ-system and for the low temperatures 
considered here. We find a mean distance between dislocations very close to 1/e 
which directly reflects the relative periodicity of the two lattices involved. 

For smaller misfits, one observes a more interesting evolution of the surface: 
initially, a pseudomorphic adsorbate film grows with the same mean lateral atomic 
distances aj at oc <j s as in the substrate. At a rather well defined film thickness, 
misfit dislocations appear everywhere on the surface, which are later overgrown. 
Snapshots of the surface evolution in the vicinity of a single dislocation are shown 
in Figure for e = —5%. 

In our most recent studies of the low misfit regime, we have slightly amended 
the model in comparison with Section [5] Deposited particles perform a so-called 
downhill funneling upon arrival at the surface: It is assumed that its kinetic energy 
enables an arriving atom to move to the lowest position in the vicinity of the depo- 
sition site. Similar incorporation effects, which are clearly not thermally activated, 
are discussed in jS] (this volume) and in greater detail in ^ [2] . The funneling 
process smoothens the surface by reducing the effect of the strong Schwoebel bar- 
rier in the system. It hence allows for the simulation of thicker films without a 
pronounced mound structure as, for instance, in Fig. [21 Note, however, that our 
findings are to a large extent robust with respect to such modifications. 



8 



Biehl, Much, and Vey 




Figure 3. The emergence of a misfit dislocation in the model 
without downhill funneling for e = —5%, other details as 
in Fig. [5J The snapshots correspond to a mean film thickness of 
12ML, 13ML and 18ML, respectively. 

For positive misfits the adsorbate is laterally compressed and it reacts by 
assuming a vertical lattice spacing a± which is larger than in its undisturbed 
bulk. A simple consideration yields for small e in the LJ-system a vertical spacing 
^const a; ot -y/3/4(l + 4/3e) of the atomic layers in the pseudomorphic film. After 
misfit dislocations have emerged, the adsorbate approaches its relaxed undisturbed 
bulk structure, as indicated by the decrease of a± with the film thickness. Figure^] 
shows the result of our simulations on average over 10 independent runs for system 
sizes of at least L = 400 and misfits in the range 1.4% < e < 2.2%. The qualitative 
evolution of a± agrees with novel measurements of this quantity in experimental 
studies of II-VI-semiconductor heteroepitaxy |29| . In a forthcoming publication 
the comparison with experimental data will be presented in greater detail. 

Perhaps the most interesting result in this context concerns the critical film 
thickness h c at which the dislocations appear in the system. In Fig.QJ it is marked 
by the deviation of a± from its initial value. The rescaling of the thickness with 
e -3/2 m Fig. ^ shows a relatively good collapse of the curves in the relevant phase 
of growth. It is by no means our intention to suggest that, here, a true dynamical 
scaling law exists as in, e.g., kinetic roughening The small range of consid- 

ered misfits would certainly not allow for such a claim. Nevertheless, our data is 
consistent with a critical thickness of the form 

h c cx e~ 3 / 2 . (3.1) 

In the previous study of a slightly modified model we considered a different measure 
of the critical film thickness JHj ■ There, we found the same power law behavior for 
positive and negative misfits in a much wider range of misfits. Preliminary studies 
yield the same qualitative result for the model with other pair potentials in place 
of the Lennard- Jones interactions. Hence we believe that the observed dependence 
is robust with respect to details of the model and displays a certain degree of 
universality. 
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Figure 4. Dislocations: development of the vertical layer spac- 
ing in the model with downhill tunneling, see text, with increasing 
film height for various values of the misfit. Results were obtained 
on average over 10 independent simulation runs and the spacing 
Si is rescaled by its initial value. The scaling of the film thickness 
is according to Eq. (|3.1|) and results in a relatively good collapse 
of the curves close to the critical thickness. 

Note that the power law behavior l|3.1fl disagrees with the results of energy or 
force balance considerations, as for instance the well-known and widely accepted 
relation derived by Matthews and Blakeslee [SUj. One can argue, however, that the 
dislocation formation as observed here is a kinetic phenomenon far from equilib- 
rium. Indeed, alternative approaches as suggested by Cohen-Solal and co-workers 
predict the power law dependence for the critical thickness pH] EH] ■ Further- 
more, experimental data for several IV-IV, III-V, and II- VI semiconductor systems 
shows very good agreement with the hypothesized power law 

In forthcoming investigations we will study the dependence of the surface 
evolution on the temperature and deposition rate. This should provide further 
evidence for the robustness of relation (|3.1() and for the interpretation of misfit 
induced dislocation formation as an activated process. 

4. Stranski-Krastanov like growth 

Obviously, dislocations should dominate the strain relaxation in sufficiently thick 
films and for large misfits. In material systems with relatively small mismatch an 
alternative effect governs the initial growth of very thin films: Instead of growing 
layer by layer, the adsorbate aggregates in three-dimensional structures, allowing 
for partial relaxation. The term 3D-islands is used to indicate that these struc- 
tures are spatially separated. The effect is to be distinguished from the emergence 
of mounds due to the Ehrlich-Schwoebel instability ^ |2] which also occurs in 
homoepitaxy. 
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At least two distinct growth scenarios display 3D-island formation: In Volmer- 
Weber growth, the structures emerge immediately upon the substrate before the 
adsorbate even forms a closed layer P^. The situation resembles the formation of 
non-wetting droplets of liquid on a surface. It is often observed in systems where 
adsorbate and substrate are fundamentally different, an example being Pb on a 
graphite substrate 

Here, we will focus on on the so-called Stranski-Krastanov (SK) growth 
mode, where 3D-islands are found upon a persistent pseudomorphic wetting-layer 
(WL) of adsorbate material ^[S]. Most prominent examples for SK-systems are 
Ge/Si and InAs/GaAs where, as in almost all cases discussed in the literature, the 
adsorbate is under compression in the WL. 

In order to avoid conflicts with more detailed or more restricted definitions of 
SK growth in the literature, we use the term SK-like growth here. It summarizes 
the following sequence of events during the deposition of a few monolayers (ML) 
of material: 

1. Initial layer by layer growth of a pseudomorphic, compressed adsorbate WL. 

2. The sudden appearance of 3D-islands, marking the so-called 2D-3D- or SK- 
transition at a kinetic WL thickness ft^, L . 

3. Further growth of the 3D-islands, which is fed by additional deposition and 
by incorporation of surrounding WL atoms. 

4. The observation of separated 3D-islands of similar shapes and sizes, on top 
of a WL with reduced stationary thickness /iwl < ^wl- 

In order to avoid confusion, and since one might interpret our (l+l)-dim. 
model as a cross section of the full (2+l)-dim. picture, we will use the term 2D- 
3D-transition as usual throughout the following. 

In SK growth a number of effects might play important roles, including the 
mixing interdiffusion of materials or the segregation of compound adsorbates. 
These effects are certainly highly relevant in many cases, see several contribu- 
tions in 0] and 0. However, SK-like growth is observed in a variety of material 
systems which may or may not display these specific features. For instance, inter- 
mixing or segregation should be irrelevant in the somewhat exotic case of large 
organic molecules like PTCDA deposited on a metal substrate. Nevertheless, this 
system shows SK-like growth according to the above definition • 

This very diversity of SK-systems gives rise to the hope that this growth 
scenario might be explained in terms of a few basic mechanisms. Accordingly, it 
should be possible to capture and identify these universal features in relatively 
simple model systems without aiming at the reproduction of material specific 
details. This hope motivated the investigation of SK-growth in the frame of our 
off-lattice model. 

An important modification beyond the description in Section [2] concerns the 
interlayer diffusion of adatoms. As argued above, the Ehrlich-Schwoebel effect 
would be much less pronounced in the physical (2 + l)-dim. situation. In our 
investigation of the SK-like scenario we remove the ES-barrier for all interlayer 
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Figure 5. SK-like growth: A section of a simulated crystal 
with six layers of substrate as obtained for e = +4% and 
Rd = 7.0ML/s at T = 50QK after deposition of (top to bottom) 
1.5ML, 2.3ML, 3.QML and 4.0ML. In the second snapshot the 
critical WL thickness h\ VL sw 2.3ML is reached and island forma- 
tion sets in. Note how smaller islands grow and larger structures 
break up. Note also that WL particles contribute to the growth 
of islands. Eventually, the well separated structures are located 
on a stationary WL with, here, hyyL ~ 1- The darker a particle 
is displayed, the larger is the average distance from its nearest 
neighbors. White represents a s whereas a a corresponds to a dark 
grey on this scale. 

diffusion events at terrace edges by hand. One motivation is the above mentioned 
over-estimation. More importantly, we wish to investigate strain induced island 
formation without interference of the ES instability. Note that the latter leads to 
the formation of mounds even in homoepitaxy ^ [2] . 

In order to favor the emergence of a wetting layer, we set U s > U as > 
U a in our model. As an example we have used U s = l.OeV^, U a = 0.74eV, and 
U as = Q-86eV accordingly. If not otherwise specified, we consider a positive misfit 
of e = 4% in the following. We have demonstrated before J2| j cf. Section 
that strain relaxation through dislocations is not expected for e = 4% within 
the first few adsorbate layers. Indeed, no misfit dislocations were observed in the 
simulations presented here. 

In the simulation we realize a situation very similar to many experiments: 
a fixed amount of adsorbate material, corresponding to 4ML, is deposited at a 
constant rate Rd- After deposition ends, we allow for a short relaxation period, 
in which atoms can still diffuse on the surface. Note however, that the essential 
surface properties are already determined during growth. 
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Our model displays a behavior along the lines of our operative definition of 
SK-like growth. The snapshots of an example simulation run in FigureJSJcorrespond 
precisely to the four stages outlined above, illustrating mpeg-movies of sample 
simulations are available at |35| . 

Several properties of surface diffusion in our model are discussed in |21j in 
greater detail. Here, we shortly summarize the essential features: 

a) Adsorbate diffusion right on the substrate is relatively slow, due to the fact 
that U as is quite large and satisfies U as > U a . 

b) Adsorbate diffusion on a strained WL is relatively fast. The diffusion barrier 
decreases with the thickness, but essentially saturates at 3 or AML in our 
model, where the interaction with the substrate becomes negligible. 

c) An adsorbate particle on top of an existing island is subject to a strong 
diffusion bias towards the island center. This result was already reported in 
|lfi| . The bias is due to the partial relaxation in the island and unrelated to 
the Ehrlich-Schwoebel effect, which has been eliminated in our model. 

The first two effects, (a) and (b), clearly favor and stabilize the existence of 
a wetting layer. Qualitatively the same relation is obtained in experimental inves- 
tigations of the Ge/Si system an d 0. With the model parameters specified 
above we find an activation barrier of approximately 0.57eV for directly on the 
substrate and roughly 0A7eV for diffusion on the first wetting layer. 

On the contrary, the diffusion bias (c) stabilizes existing islands by essentially 
confining adatoms to their top terraces. Of course this is the case for particles de- 
posited onto the island. More importantly, we find a significant contribution of 
particles that jump upward, from the WL atop an island. The rate for such pro- 
cesses is quite small, generally, but becomes significant close to the SK-transition, 
see [21] for details. 

In simulations with different deposition rates we observe that the kinetic WL 
thickness increases with Rd [213 EH If the formation of second or third layer nuclei 
by freshly deposited particles was the driving force, one would expect more frequent 
nucleation and an earlier 2D-3D-transition at higher growth rates. We conclude 
that, in the main, upward jumps trigger the SK-transition. Further evidence and 
a more detailed discussion of this important point is given in |21j . 

Our investigations suggest the following picture of the Stranski-Krastanov 
transition: The diffusion properties, (a) and (b), favor the formation of a wetting 
layer and stabilize it. As the film grows, strain accumulates in the film and upward 
hops from the WL become more probable. The precise spatial modulation of Rd in 
the film presumably determines the final island sizes and will be studied in forth- 
coming investigations. Once multilayer islands have emerged they are stabilized 
by the diffusion bias (c). 

After the 2D-3D-transition, islands grow by incorporating newly deposited 
material, but also by consumption of the surrounding WL. Note that very large 
islands are observed to split by means of upward diffusion events onto their top 
layer, cf. Fig. [SJ The migration of WL particles towards and onto the islands 
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Figure 6. Average base size (b) of multilayer islands as a func- 
tion of Rd at T — 480 and T = 500K, together with standard 
error bars. Results were obtained on average over 15 independent 
simulation runs. The inset shows the result for T = BOOK, Rd = 
A.hML/s and different misfit parameters, the solid line corre- 
sponds to (b) = 0.91/e. 



can well extend into the relaxation period after deposition ends. Eventually, a 
stationary WL thickness /iwl ~ 1 is observed in our example scenario with U as ~ 
0.86eV. By increasing the strength of the adsorbate/substrate interaction we can 
achieve, e.g., hwh ~ 2 for U as » 2.76^, but it is difficult to stabilize a greater 
stationary thickness. This effect is related to the effective short range nature of the 
LJ-potential, which is very weak for distances larger than 3tr. With, e.g., long-range 
or multi-particle interactions, the model should yield greater WL thicknesses. 

Finally, we discuss some properties of the emerging islands or SK-Dots. Fig- 
ure El displays the average lateral size, measured as the number of particles in the 
island bottom layer. The results shown here were obtained at the end of a short re- 
laxation period with Rd — 0. Whereas the mean values do not change significantly, 
fluctuations are observed to decrease in this phase. 

We observe for two different temperatures that the island size decreases with 
increasing deposition rate. This is in accordance with several experimental obser- 
vations [SZj. However, the size becomes constant and independent of T for large 
enough deposition flux. A corresponding behavior is found for the island density 
and their lateral spacing, which hints at a considerable degree of spatial ordering 
[201 . This saturation behavior further demonstrates the relevance of upward 
hops. Standard arguments |2] show that a dominant aggregation of deposited 
particles on top terraces would yield an island density that continues to increase 
with Rd- 

The inset in Fig. displays the mean island size for T = 500fsT and Rd = 
A.bML/s, i.e. in the saturation regime. Our result is consistent with a simple power 
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law of the form (b) oc 1/e. Very far from equilibrium, the only relevant length scale 
in the system appears to be, again, the relative periodicity 1/e of adsorbate and 
substrate lattice. 

5. Summary and outlook 

Despite its conceptual simplicity and the small number of parameters, our model 
reproduces several features of heteroepitaxial growth. Strain effects are incorpo- 
rated in a natural fashion as they emerge directly from the interaction of particles. 
For instance, we have demonstrated that misfit dislocations appear in the system 
when the film thickness exceeds a characteristic value. This characteristic height 
displays a power law behavior on the lattice misfit. 

We furthermore believe that, with appropriately chosen interaction param- 
eters, our model is capable of reproducing the three essential growth modes: ex- 
tended layer by layer growth (for very small misfits), Volmer- Weber for U as < U a 
and, as demonstrated already, SK-growth for U as > U a . The small number of pa- 
rameters should allow for determining the corresponding phase diagram of growth 
modes in the frame of our model system. 

Besides the exploration of the available parameter space, we intend to extend 
the model conceptually in several directions. As just one example, intermixing 
and seggregation should be included to study the relevance of these effects in the 
Quantum Dot formation. 

In order to test the potential universality of our findings, we will introduce 
different types of interactions in our model. Preliminary results for, e.g., Morse 
potentials |25| and modifications of LJ-potentials confirm our results in the context 
of dislocation formation. Ultimately, we will extend our model to the physical case 
of growth in (2 + 1) dimensions and to realistic empirical potentials for metals or 
semi-conductors. As a first step, we are currently investigating the sub-monolayer 
regime in strained heteroepitaxy of fee materials. 

In a sense, the off-lattice KMC method provides a link between traditional, 
lattice based methods and Molecular Dynamics. Hence, it might play a significant 
role in the further development of the multiscale approach. 
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